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16.  Abstract 

A  mathematical  model  is  being  developed  that  will  permit  predictions  of  the  strength  of  fiber 
reinforced  composites  containing  known  flaws  to  be  made  from  the  basic  properties  of  their 
constituents.  The  approach  is  to  embed  a  local  heterogeneous  region  (LHR)  surrounding  the  crack 
tip  into  an  anisotropic  elastic  continuum.  The  intent  is  to  have  the  model  (1)  permit  an 
explicit  analysis  of  the  micromechanical  processes  involved  in  the  fracture  process,  and 
(2)  remain  simple  enough  to  be  useful  in  practical  computations.  Material  properties  used  in 
the  analysis  are  to  be  obtained  from  a  concurrent  experimental  program  being  carried  out  at 
Virginia  Polytechnic  Institute  and  State  University. 

Computations  for  arbitrary  flaw  size  and  orientation  under  arbitrary  applied  load  combinations 
have  been  performed  for  unidirectional  composites  with  linear  elastic-brittle  constituent 
behavior.  The  mechanical  properties  were  nominally  those  of  graphite  epoxy.  With  the  rupture 
properties  arbitrarily  varied  to  test  the  capability  of  the  model  to  reflect  real  fracture  modes 
in  fiber  composites,  it  is  shown  in  the  report  that  fiber  breakage,  matrix  crazing,  crack 
bridging,  matrix-fiber  debonding,  and  axial  splitting  can  all  occur  during  a  period  of  (gradually) 
increasing  load  prior  to  catastrophic  fracture.  Of  most  importance,  the  computations  reveal 
qualitatively  the  sequential  nature  of  the  stable  crack  growth  process  that  precedes  fracture 
in  composites.  Quantitative  comparisons  with  the  VPISU  experimental  results  on  edge-notched 
unidirectional  graphite  epoxy  specimens  have  also  been  made  and  were  found  to  be  in  fair 
agreement. 
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FUNDAMENTAL  ANALYSIS  OF  THE  FAILURE  OF 
POLYMER-BASED  FIBER  REINFORCED  COMPOSITES 


INTRODUCTION 


There  are  many  different  ways  in  which  a  structure  made  of  a  fiber 
reinforced  composite  material  can  become  unable  to  adequately  perform  its 
primary  function.  In  each  such  instance  failure  is  considered  to  have  oc¬ 
curred.  The  possible  failure  modes  therefore  encompass  a  wide  range  of  pos¬ 
sibilities  from  simple  loss  of  structural  rigidity  due  to  gross  inelastic  de¬ 
formation  (e.g.,  yielding),  through  a  reduction  in  load-carrying  capacity 
due  to  localized  damage  or  separation  (e.g.,  interply  delamination) ,  to  the 
complete  loss  of  strength  due  to  large-scale  crack  growth  and  fracture.  Each 
of  these  failure  modes  can  be  gradual  or  rapid  depending  on  the  nature  of  the 
applied  loads,  the  material  properties,  the  geometry  of  the  structure,  and  the 
presence  of  cracks  or  flaws.  For  polymer-based  composites,  the  loading  rate, 
the  temperature,  and  previous  load  history  can  also  play  prominent  roles. 

In  the  work  described  in  this  report,  emphasis  is  placed  upon 
fracture  and,  therefore,  the  work  will  be  primarily  concerned  with  the  "strength” 
of  fiber  composites  containing  known  flaws.  The  term  strength  is  convention¬ 
ally  taken  to  mean  the  load  level  at  which  failure  occurs  in  a  standard  test 
specimen.  Clearly,  the  strength  will  be  a  function  of  many  different  para 
meters  arising  in  the  test  program  and  may  or  may  not  be  directly  applicable 
to  engineering  design  situations.  The  primary  purpose  of  the  research  des¬ 
cribed  in  this  report  is  to  provide  a  bridge  between  standard  laboratory  test 
procedures  and  actual  engineering  applications  of  fiber  composites  that  will 
allow  accurate  reliable  estimates  of  the  failure  loads  for  aircraft  and  other 
engineering  structures.  |  Such  a  capability  does  not  presently  exist. 

The  design  tool  that  should  be  developed  for  the  safe  and  efficient 
utilization  of  composite  materials  is  a  predictive  capability  for  fracture 
that  can  take  account  of  the  applied  loading,  the  geometry  of  the  structure, 
and  the  environmental  effects  in  terms  of  readily  measurable  properties  of  the 
composite ^s  constituents  and  of  its  microstructural  design.  The  primary 
benefits  accruing  from  such  an  analytical  capability  are  twofold.  First, 
for  a  given  structure,  material,  and  load,  the  critical  flaw  size  and  type  in 
a  structural  component  can  be  estimated  for  comparison  with  actual  flaws  ob¬ 
served  in  an  inspection  program.  Second,  guidelines  can  be  provided  tor  tne 
designer  to  tailor  more  fracture  resistant  "tougher”  materials  for  specific 
engineering  applications. 
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PROGRAM  PLAN  AND  SULDIARY  OF  PROGRESS 


The  work  described  in  this  report  represents  the  first  year  of  effort 
in  what  is  planned  as  a  multiyear  program.  The  specific  objective  of  the 
first  year’s  work  was  to  develop  the  basic  mathematical  procedures  required 
for  the  analysis  model.  To  implement  this  work,  attention  was  focussed  on 
unidirectional  composites  with  linear  elastic-brittle  material  behavior.  In 
subsequent  work,  the  model  will  be  extended  to  treat  angle  ply  laminates  and 
will  include  further  refinements  (e.g.,  inelastic  constituent  behavior)  re¬ 
quired  in  order  to  treat  actual  engineering  problems. 

The  primary  objective  of  the  work  will  be  achieved  only  if  the 
mathematical  model  developed  is  capable  of  delineating  the  role  of  the  various 
micromechanical  failure  processes  that  dictate  the  ultimate  failure  point  of 
fiber  reinforced  composites.  The  research  described  in  this  report  seeks  this 
end  by  merging  a  micromechanical  failure  analysis  with  a  macromechanical  frac¬ 
ture  mechanics  approach.  This  approach  treats  the  material  as  heterogeneous 
and  anisotropic  where  microstructural  effects  predominate  and  as  homogeneous 
and  anisotropic  where  it  is  permissible  and  practical  to  do  so.  In  this  way, 
direct  consideration  can  be  given  to 

•  The  external  size  and  shape  of  the  structure  and  the 
laminate  stacking  sequence 

•  The  applied  loads  acting  on  the  structure,  both 
mechanical  and  thermal,  and  environmental  effects 

•  The  size,  shape,  and  orientation  of  a  flaw  in  the 
laminate. 

In  particular,  the  manner  in  which  these  parameters  influence  the  sequence 
of  microstructural  failure  events  whereby  a  flaw  extends  stably  under  an 
increasing  load  up  to  the  point  of  catastrophic  fracture  will  be  determined. 

It  should  be  understood  that  the  conventional  approach  to  fracture  analysis, 
linear  elastic  fracture  mechanics,  is  not  capable  of  coping  with  this  degree 
of  complexity,  thus  necessitating  the  more  general  development  being  pursued 


here . 
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Linear  elastic  fracture  mechanics  (LEFM)  is  a  predictive  technique 
applicable  to  structural  components  containing  crack-like  flaws  when  the 
material  used  fits  certain  key  assumptions  used  in  the  LEFM  theory.  Specifi¬ 
cally,  the  material  must  behave  very  nearly  as  would  a  completely  linear 
elastic  perfectly  homogeneous  ideal  material.  The  LEFM  approach  has  been 
successful  when  properly  applied  (e.g.,  to  high  strength  steel),  but  consider¬ 
ably  less  successful  in  applications  where  the  basic  assumptions  are  not  well 
satisfied.  Fiber  reinforced  composite  materials  are  an  outstanding  example 
of  the  latter  case.  In  particular,  in  fiber  composites  with  a  flaw,  a  non¬ 
linear  "damage”  zone  is  generally  produced  at  the  flaw.  This  zone,  generally 
growing  in  a  stable  manner  under  an  increasing  applied  load,  has  a  profound 
effect  on  the  eventual  point  of  complete  fracture.  This  is  shown,  for  example, 
in  the  work  of  Mandell,  et  al.  [l]*  as  well  as  by  many  other  investigators « 

The  damage  zones  in  a  fiber  composite  are  the  result  of  a  large 
number  of  discrete  failure  processes,  e.g.,  fiber  breaking,  matrix  yielding, 
etc. ,  that  occur  in  the  highly  stressed  region  ahead  of  a  crack-like  flaw. 

These  individual  micromechanical  events  do  not  conform  to  the  basic  assumptions 
of  LEFM.  Consequently,  it  is  not  surprising  that  the  agregate  of  such  pro¬ 
cesses  cannot  be  treated  by  LEFM.  What  is  therefore  needed  is  a  generalized 
fracture  mechanics  treatment  which  explicitly  recognizes  the  fundamental  na¬ 
ture  of  damage  growth  in  composite  materials.  Specifically,  a  proper  analysis 
of  fracture  in  fiber  composites  must  be  cognizant  of  two  key  features:  the 
generally  anisotropic  constitutive  behavior  of  the  material,  and  the  hetero¬ 
geneous  nature  of  the  fracture  process. 

The  approach  developed  in  this  report  can  be  likened  conceptually 
to  boundary  layer  theory  and,  in  application,  to  the  well-established  singu¬ 
lar  perturbation  and  matched  asymptotic  expansion  techniques  of  fluid 
mechanics.  That  is,  the  problem  of  a  composite  material  containing  a  flaw 
is  divided  into  distinct  "inner"  and  "outer"  regions.  In  each  of  these  regions, 
the  material  is  modeled  in  different  ways.  The  inner  region,  which  contains 
the  crack  tip,  is  considered  on  the  microscopic  level  and  treats  the  material 
as  being  heterogeneous.  The  outer  region  surrounds  the  crack  tip  region.  It 
is  treated  as  a  homogeneous  orthotropic  continuum.  For  simplicity,  the  inner 
crack  tip  region  will  be  referred  to  in  this  report  as  the  LHR  (for  local  hetero 
geneous  region)  while  the  outer  region  will  be  simply  called  the  continuum. 

*  References  are  given  on  page  ''^5. 
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The  LHR  consists  of  elements  representing  the  matrix,  the  fibers, 
and  the  fiber-matrix  interface  zones.  Ttie  constitutive  relations  of  these 
elements,  up  to  and  including  their  rupture  points,  are  presumed  to  be  knoxvrn 
from  experiments.  Any  element  of  a  fiber  composite  ruptures  when  an  intrinsic 
critical  energy  dissipation  rate  can  be  provided  at  some  point  of  that  element. 
These  critical  values  are  assumed  to  be  independent  of  the  local  stress  field 
environment  at  the  point  of  incipient  rupture.  This  permits  data  from  fracture 
tests  on  isolated  fibers,  matrix  material,  and,  possibly,  unflawed  composites 
(to  obtain  interface  strengths)  to  be  directly  inserted  into  the  model.  Ma¬ 
terial  properties  used  in  the  analysis  work  and,  ultimately,  critical  tests 
of  the  predictions  of  the  model  will  be  obtained  from  a  concurrent  experi¬ 
mental  program  being  carried  out  under  NASA  sponsorship  at  Virginia  Polytechnic 
Institute  and  State  University. 

Progress  to  date  has  permitted  computations  to  be  performed  for 
unidirectional  composites  with  e las t ic- per feet ly  brittle  constituent 
behavior.  The  mechanical  properties  have  been  those  of  graphite  epoxy.  The 
rupture  properties  have  also  been  arbitrarily  varied  to  test  the  capability 
of  the  model  to  reflect  real  fracture  modes  in  fiber  composites.  It  has 
been  shown  that  fiber  breakage,  crack  bridging,  matrix-fiber  debonding,  and 
axial  splitting  can  all  occur  during  a  period  of  (gradually)  increasing  load 
prior  to  catastrophic  fracture.  In  this  way,  the  sequential  and  interrelated 
manner  in  which  each  individual  local  rupture  event  occurs  during  damage 
growth  preceding  fracture  in  fiber  composites  is  revealed  by  the  computations. 
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DESCRIPTION  OF  THE  ANALYSIS  TECHNIQUE 

In  the  preceding  section  of  this  report,  the  objectives  of  the  work 
and  an  outline  of  the  general  approach  were  given.  The  focal  point  for  this 
description  was  the  two-dimensional  local  heterogeneous  region  (LHR)  surrounding 
the  crack  tip.  A  typical  LHR  model  is  shown  in  Figure  1.  Depicted  is  a  uni¬ 
directional  fiber  composite  containing  three  distinct  components:  the  fibers, 
the  matrix,  and  the  fiber-matrix  interface  zones.  In  this  section  of  the  report, 
the  discrete  elements  that  comprise  the  LHR  and  the  manner  in  which  they  are 
assembled  to  give  a  quantitative  predictive  capability  for  composite  fracture 
will  be  described  in  detail. 


LHR  Boundary  Conditions 


In  the  work  performed  so  far,  the  interaction  between  the  LHR  and  the 
continuum,  as  is  appropriate  for  preliminary  stages  of  the  work,  has  been  taken 
in  the  simplest  possible  manner.  This  is  by  specifying  the  displacements  on  the 
boundary  of  the  LHR  in  accord  with  the  "rigid  boundary  conditions"  approach.  Use 
of  this  scheme  is  tantamount  to  assuming  that  the  LHR  is  large  enough  that  the 
nonlinear  inhomogeneity  of  the  crack  tip  region  does  not  affect  the  periphery 
of  this  region.  In  other  words,  the  LHR  and  the  continuum  are  assumed  to  be 
uncoupled.  Consequently,  the  displacements  at  the  LHR  boundary  are  exactly  the 
same  as  if  the  entire  cracked  body  is  an  elastic  continuum.  The  required  rela¬ 
tions  can  therefore  be  obtained  from  the  work  of  Sih  and  Liebowitz  [2].  Their 
approach  is  briefly  summarized,  as  follows. 

Consider  a  plate  with  its  major  dimensions  lying  in  the  xy  plane. 

For  either  plane  stress  or  plane  strain  deformation,  the  inplane  strains 

£  and  Y  will  depend  only  on  the  inplane  stresses  a  ,  a  ,  and  t  . 
y  ’  xy  X  y  xy  ^ 

Hence,  the  appropriate  constitutive  relations  for  a  rectilinearly  anisotropic 
body  in  a  state  of  plane  deformation  are  given  by 


=  . 

a 

+ 

a 

a, . 

T 

11 

X 

12 

y 

16 

xy 

=  CL.  r. 

a 

+ 

0 

+ 

T 

12 

X 

22 

y 

26 

xy 

=  Ci^  ^ 

a 

+ 

0 

+ 

a, . 

T 

y  16 

X 

26 

y 

66 

xy 
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Continuum 


Material  Orientation 
(Arbitrary) 
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Sih  and  Liebowitz  show  that  the  solution  to  the  governing  differential  equation 
of  two-dimensional  anisotropic  elasticity  theory  is  associated  with  the  roots 
of  the  characteristic  equation 


^11 


3  2 


2  aze  ^  +  <Xz2 


It  can  be  shown  that  the  roots  of  Equation  (2)  are  either  complex  or  are  pure 
imaginary.  Hence,  these  can  be  labeled  2S2>  3.nd  This  suggests  the 

introduction  of  the  complex  variables  =  x  +  .6^y  and  =  x  +  ^2^*  plane 

problem  of  an  anisotropic  body  is  thereby  reduced  to  the  determination  of  two 
complex  potential  functions  of  a  complex  variable  (^(z^)  and  ^(^2)  that  satisfy 
the  prescribed  boundary  conditions  of  the  problem. 

As  further  shown  in  reference  2,  the  displacement  field  is  expressed 
in  terms  of  the  potential  functions  by  the  relations 


Omitting  the  details,  potential  functions  for  a  cracked  body  infinite  in  extent 
have  been  determined  for  insertion  into  Equations  (3).  In  particular,  a  solution 
for  remote  loading  consisting  of  a  uniform  tensile  loading  a  acting  in  the  di¬ 
rection  normal  to  the  crack  plane  and  a  shear  loading  t  parallel  to  the  crack 
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plane  for  a  crack  of  length  2a  can  be  obtained.  For  a  polar  coordinate  system 
with  origin  at  the  crack  tip  as  shown  in  Figure  2,  the  displacements  near  the 
crack  tip  are  given  by 


u  =  K  (— I  Rc  I  ^  i^p^Ccos.})  +  ^2  sincji)  + 

lU/  ^2-^ 


p^  (cos  cj)  4-  sin  ) 


Re  I ^2  <t>  +  ^2 


I/2I  f 

-  ^2  ^2^  (cos  (p  +  6^  sin  b  )  V  + 


'’2  <!>  +  ^2  + 


(cos  ^  ^  ) 
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where 


-  o  /'tt  a 


and 


(6) 


=  T  /tt  a 

are  the  Mode  I  and  Mode  II  stress-intensity  factors  for  the  problem. 

Equations  (4)  and  (5)  can  be  used  to  determine  displacements  on  the 
LHR  boundary  via  the  rigid  boundary-condition  approach.  This  means  that  the 
applied  stresses  acting  on  the  body  are  transmitted  through  the  continuum  region 
to  the  crack-tip  region  and  are  "sensed”  at  the  crack  tip  in  terms  of  the  stress- 
intensity  factors  and  An  independent  specification  of  the  load  and 

crack  length  is  therefore  unnecessary.  Hence,  although  derived  for  an  infinite 
medium,  the  approach  can  be  used  for  bodies  with  finite  boundaries  by  simply  in¬ 
serting  the  appropriate  stress-intensity  factors.  In  doing  this,  it  must  be 
tacitly  assumed  that  the  LHR  is  (1)  large  enough  relative  to  the  micros t rue tural 
dimensions  of  the  composite  that  the  boundary  displacements  are  closely  given 
by  continuum  theory  and  (2)  small  enough  relative  to  the  crack  length  and  di¬ 
mensions  of  the  body  that  the  singular  behavior  of  the  continuum  solution  at  the 
crack  tip  dominates. 

While  appropriate  for  preliminary  work,  it  is  anticipated  that  the 
rigid  boundary  condition  approach — even  with  periodic,  updating  to  reflect  the 
progress  of  the  crack  through  the  LHR — may  prove  to  be  too  restrictive.  There¬ 
for,  in  subsequent  work,  a  flexible  boundary-condition  approach  which  extends 
an  approach  developed  in  work  previously  carried  out  at  Battelle  [3]  will  be 
used.  This  is  described  in  more  detail  in  the  Recommended  Further  Research 
section  of  this  report.  It  might  be  noted  that  xvith  flexible  boundary  conditions, 
crack  length  and  load  must  be  specified  individually. 
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The  LHR  Element 


As  shown  in  Figure  1,  the  LHR  for  a  fiber  reinforced  composite  is 
considered  to  be  made  up  of  three  different  types  of  constituents:  fiber, 
matrix,  and  fiber-matrix  interface.  Each  of  the  individual  constituents  in 
the  LHR  must  be  capable  of  rupturing  to  allow  the  body  to  exhibit  the  changes 
in  strength  that  correspond  to  various  levels  and  orientations  of  local  break¬ 
age.  A  good  deal  of  success  has  been  obtained  by  Kanninen  [4,5]  using  "spring¬ 
like"  elements  to  model  a  local  fracture  phenomena.  This  fact,  taken  together 
with  the  increased  computational  simplicity  of  this  approach,  has  led  to  the 
adoption  of  the  basic  element  shown  in  Figure  3  for  this  work. 

As  shown  in  Figure  3,  each  element  has  eight  degrees  of  freedom 
and  is  connected  in  the  LHR  at  its  four  corner  node  points.  Extensional  stiff¬ 
ness  is  provided  by  four  extensional  connectors.  These  connectors  resemble 
simple  extensional  springs,  but  also  have  a  lateral  contraction  or  Poisson  effect. 
A  material  made  up  of  a  set  of  these  connectors  behaves  in  uniform  extension 
exactly  like  a  homogeneous  orthotropic  material.  Similarly,  shear  stiffness 
is  incorporated  in  the  element  through  a  rotary  spring  at  each  node  point  to 
give  the  proper  response  to  shear  loadings.  The  values  of  the  spring  con¬ 
stants  are  functions  of  the  material’s  elastic  constants  and  of  the  element 
size  and  shape. 

In  addition  to  giving  the  proper  response  to  loads,  the  LHR  elements 
fracture  according  to  the  "codes"  shown  in  Figure  4.  In  Figure  4,  Fracture  Code 
1  represents  an  increment  of  crack  extension  in  the  plane  lying  midway  between 
Nodes  1  and  4  that  extends  from  the  left  side  of  the  element  to  the  center  of 
the  element.  Fracture  Code  *2  represents  the  continuation  of  the  crack  to  the 
right  side  of  the  element.  Fracture  Codes  3  and  4  represent  similar  increments 
of  cracking  in  the  vertical  direction. 

Using  an  energy  approach,  it  is  possible  to  trace  those  components 
of  the  elemental  stiffness  matrix  that  are  associated  with  each  fracture  code  in 
every  element  in  the  LHR.  This  information  can  be  assembled  into  the  stiffness 
matrix  such  that  a  solution  can  be  obtained  for  the  incipient  rupture  of  any 
subelement  in  the  LHR.  Note  that  this  is  most  easily  possible  if  (as  in  the 
present  work),  the  constituent  behavior  is  restricted  to  being  linear  elastic 
to  the  rupture  point.  Inelastic  behavior,  while  more  complicated,  is  not 
precluded,  however. 
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Knowledge  of  the  stiffness  components  attributable  to  each  fracture 
condition  also  provides  a  method  of  calculating  the  energy-release  rate  for 
crack  advance  by  any  of  the  four  codes  provided  for  in  each  LHR  element.  Of 
course,  the  critical  rupture  energy  values  must  be  specified  to  provide  a  de¬ 
cision  rule  for  breakage  in  each  separate  element.  Note  that  while  the  model 
allows  separate  critical  values  to  be  specified  for  each  of  the  four  fracture 
codes  for  every  one  of  the  elements  in  the  LHPv,  ordinarily  different  values 
will  be  specified  only  for  the  different  constituents,  i.e,,  fiber,  matrix,  or 
interface . 

Elemental  Stiffness  Formulation 


Energy  principles  can  conveniently  be  used  to  determine  the  stiffness 
of  the  ERR  structure.  The  sign  conventions  for  an  LIIR  element  will  be  taken  as 
shoTvn  in  Figure  5.  Then,  using  the  following  notation. 


C. 


xy 


j-  l-T 

rotational  spring  constant  at  the  i-~  node 


extensional  stiffness  in  the  x  direction  betv*7een  the 

.th  tll 

i —  and  j —  nodes 

extensional  s  tit  mess  in  the  y  direction  betw'een  the 

,  til  ,  th 

1 —  and  J-—  noaes 

cross  extensional  component  (effect  of  Poisson's  ratio) 

V  ,  . ,  .  til  ,  .  th 

beeween  cne  i —  ana  i —  nodes 


u^  ~  X  displacement  at  the  i*”^  node 


V.  -  y  displacement  at  the  i-^  node 


the  to  La.!,  strain  energy  U  stored  in  the  element  for  any  set  of  arbitrary  nodal 

displacements  u.  and  v,  can  be  written  as 
1  1 
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This  is  accomplished  by  relating  the  node  forces  to  the  partial  derivatives  of 
the  total  strain  energy.  That  is 


and 


= 

y .  3v . 

^  1  1 


(9) 


For  example, 


""l  ^2j 

2(u  -u  )“ 

4- 

Ax  Ayl 

+  2 
(Ay)  J 

'  2(v^-V3)  2(u^-u^) 

Ax  Ay  /,  n2 

(Ay) 


(10) 


Expansion  of  this  relation  will  give  the  elements  in  the  first  row  of  the  matrix 
[K  ]  .  Subsequent  derivatives  taken  with  respect  to  the  other  degrees  of  freedom 
generate  the  remaining  rows. 

The  spring  constants  used  in  the  formulation  must  be  related  to  the 
constituent's  elastic  properties.  This  is  accomplished  via  the  relations 


‘xy 

^12 

for 

i  J 

=  1,...4 

■ij 

=  Q.,  ^ 

for 

j 

=  1,....,4 

XX 

^11  Ax 

■ij 

=  Q  — 

for 

all 

i,  j  =  1, ... ,4 

yy 

22  Ay 

i 

i  Q33 

Ay  for 

all 

i  =  1, ... ,4 
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where  [Q]  is  defined  as  the  constitutive  stiffness  niatrix  in  the  usual  maniier 
as 


r- 1  r 

1 

1 

c 

X  i 

\ 

"x 

O' 

II 

o’ 

1 

e 

y  \ 

1 

y 

!  T  1  1 

! 

Y 

1 

1 

'xy 

L  J  L,. 

J 

J 

Specifying  these  relations  is  tantamount  to  assuming  that  the.  material  in  each 
element  is  completely  homogeneous.  The  desired  heterogen  it y  arises  from  the  fact 
that  the  are  different  for  the  fiber,  matrix,  and  interface  elements. 

The  Energy-Release  Rate 

The  local  rupture  criterion  used  in  this  work  is  a  generalized  version 
of  the  ordinary  energy-release  rate  (or  crack-driving  force)  quantity  of  LEFM. 

In  formulating  it,  attention  m.ust  be  placed  on  the  complete  structure.  The 
discretized  form  of  the  equations  of  equilibrium  for  a  complete  structure  can 
be  represented  by  the  matrix  equation 


where  k  is  the  structural  stiffness  matrix,  u  is  the  nodal  displacement  vector 
and  F  the  applied  nodal  load  vector.  The  energy-release  rate  G  can  be  defined 
in  a  global  sense  by  relating  it  to  the  change  in  work  done  by  the  applied  load¬ 
ing  and  the  strain  energy  during  a  virtual  crack  extension.  This  is 


da  da  (P3) 


where 


T 

W  “  u  F 


(14) 
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is  the  work  done  by  the  applied  loading, 


is  the  strain  energy  and  a  represents  the  crack  length.  (In  these  equations,  the 

superscript  T  denotes  the  transpose.) 

When  Equations  (14)  and  (15)  are  introduced  into  Equation  (13) ,  it 

is  found  that 


dM  T 

^  ~  da  ? 


1  T  d  , 

2  ^  do.  k  % 


It  can  be  seen  from  Equation  (12)  that  the  quantity  within  the  braces  of  Equa 
tion  (16)  is  a  null  vector.  Furthermore,  if  is  independent  of  the  crack  length, 
then  the  second  term  also  vanishes.  Equation  (16)  therefore  reduces  to 


r'  _  i  T  d  . 

^  2  ^  da  ^ 


Now,  it  can  be  considered  that  the  only  contribution  to  the  matrix^  [k]  is  that 
due  to  the  stiffness  matrix  K  of  the  element  containing  the  crack  tip.  If 
denotes  the  displacement  vector  of  the  nodes  of  the  element  only,  then  Equation 
(17)  can  be  used  to  obtain  an  approximation  of  G  that  can  be  written  as 


where  [K^]  ([K^])  is  the  elemental  stiffness  matrix  before  (after)  an  extension 


Aa  of  the  crack  within  the  element.  Recognizing  that 


1  T 

2  Ji  W 
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is  the  elemental  strain  energy,  then  an  alternative  form  of  Equation  (18)  is 


A  similar  development  for  the  energy-release  rate  occurs  in  the  recently  pub¬ 
lished  papers  by  Hellen,  et  al.  [6,?]. 

Program  Solution  Techniques 

To  use  the  analysis  that  has  been  developed  thus  far,  the  material 
properties  of  the  constituents  and  the  micromechanical  geometry  must  be  specified. 
(The  material  properties  are  also  used  to  determine  the  homogeneous  orthotropic 
properties  of  the  bulk  composite  through  an  effective  modulus  or  volume  fraction 
technique.)  The  way  in  which  the  constituent  properties  and  the  geometry  are 
used  to  generate  elemental  stiffnesses  for  assemblage  into  the  LHR  will  be  des¬ 
cribed  next. 

There  are  two  classes  of  nodal  points  used  in  the  LHR.  The  first 
consists  of  the  points  on  the  boundary  of  the  LHR  where  displacements  are  pre¬ 
scribed.  The  nodal  points  in  the  interior  of  the  LHR  having  unknown  displacements 
make  up  the  second  class.  It  is  useful  to  assemble  the  stiffness  matrix  in 
such  a  way  that  these  two  groups  are  isolated.  The  partition  used  here  is 


(21) 


In  Equation  (21),  the  subscripts  do  not  refer  to  nodes,  but  to  the  nodal  point 
classification  just  described.  That  is  (U2)  represents  the  prescribed  displace¬ 
ment  vector  of  the  peripheral  points  while  (u^)  is  the  vector  of  the  unknown  dis¬ 
placements  for  which  a  solution  is  desired.  Note  that  the  number  1  was  not  used 
as  a  subscript  because  it  is  more  convenient  to  reserve  it  for  specified  zero 
displacements . 

By  performing  the  matrix  multiplication  indicated  in  Equation  (21), 
the  following  two  equations  are  obtained: 
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[^22]  [’^23]  {“3}  ^^2^ 


(22) 


[^32]  {'^2}'*’  [^33]  l'^3}  I  ^3 1 


(23) 


Most  often,  it  is  either  impossible  or  inefficient  to  store  large  matrices  such 
as  [k^^]  in  the  computer’s  central  processor  core  area.  For  this  reason  [k^^] 
is  stored  externally.  The  external  storage  can  be  minimized  by  utilizing  the 
fact  that  the  [k^^]  matrix  is  always  symmetric  and  banded;  the  latter  term 
indicating  the  tendency  for  nonzero  values  to  cluster  around  the  diagonal  of 
the  matrix.  By  taking  advantage  of  matrix  partitioning  (i.e. ,  the  grouping 
of  prescribed  and  unprescribed  degrees  of  freedom),  symmetry,  and  the  banded 
nature  of  the  stiffness  matrix,  storage  of  only  a  small  portion  of  the  original 
total  stiffness  matrix  is  required.  However,  this  reduced  storage  requirement 
can  still  be  very  large  in  some  problems. 

Solution  of  Equation  (23)  is  accomplished  through  the  use  of 
Gaussian  elimination  with  back  substitution.  Small  parts  of  stored  [k^^]  matrix 
are  brought  into  core  as  required  in  this  operation.  Assemblage  of  the 

>v 

required  parts  of  the  stiffness  matrix  is  done  by  the  direct  stiffness  method. 
Components  that  would  ordinarily  fall  within  the  1^22^  ^^23^  submatrices 

are  discarded.  Stiffness  components  that  fall  within  the  [k^2^  submatrix  are 
used  in  a  nonzero  {F^}  vector.  Finally,  those  values  of  [k^^]  falling  on  or 
above  its  diagonal  are  stored  externally. 

Once  the  boundary  conditions  have  been  prescribed,  an  out-of-core 
type  solution  of  the  partitioned  LHR  stiffness  matrix  with  boundary  conditions 
for  a  unit  load  is  performed.  The  resulting  displacement  vector  is  then  used 
to  calculate  the  potential  fracture  energies  associated  with  each  of  the  four 
fracture  codes.  Ratios  of  the  prescribed  critical  energy  levels  to  these  cal¬ 
culated  energies  are  next  computed.  The  applied  load  is  then  adjusted  such 
that  the  highest  ratio  becomes  exactly  equal  to  unity.  This  critical  region 
(if  one  exists)  is  allowed  to  break  in  one  of  the  four  ways  (codes)  described 
above.  Appropriate  modifications  are  then  made  to  the  LHR  stiffness  matrix. 

*  A  detailed  description  of  the  direct  stiffness  assemblage  technique  used 
in  this  work  is  given  by  Desai  and  Abel  [8]. 


Because  the  system  is  complete3-y  linear  elastic,  a  solution  for  another  load  level 
can  then  be  performed  in  exactly  the  same  manner «  In  this  may^  the  properties 
of  a  crack  tip  damage  zona  as  a  function  of  an  increasing  applied  load  can  be 
generated. 

As  each  local  rupture  occurs,  the  prescribed  boundary  conditions  must 
be  adjusted.  If  the  crack  has  propagated  in  a  self--similar  fashion^  this  is 
readily  done  by  shifting  the  origin  (cf.  figure  2).  If  not.  this  may  present 
some  very  real  dif f iciilties  for  then  the  continuum  solution  no  longer  exactly 
applies  and  some  approximations  must  be  used  for  setting  the  boundary  conditions. 
The  larger  the  LHR.  the  less  likely  this  mill  be  necessary.  Of  course ^  large 
LHR's  require  large  solution  times  even  though  storage  is  generally  not  a  problem 
because  of  the  out-of-core  solution.  It  is  likely  that  the  boundary  condition 
problem,  can  be  handled  by  introducing  a  flexible  boundary-condition  scheme  (see 
Recomm-ended  Future  Research  section) g  or^  possibly 5  by  allowing  the  fractured  LHR 
to  interact  in  a  hybrid  continuum  model  for  future  boundary  updating. 

Verification  of  Computational  Model 

Before  performing  extensive  calculations  on  composite  materials g  it 
is  appropriate  to  verify  the  computational  model  by  checking  it  against  known 
solutions.  B}/'  letting  the  model  represent  a  linear  elastic  homogeneous  material ^ 
two  checks  can  be  made.  These  are  on  the  displacements  in  the  LHR  and  on  the 
calculated  energy-release  rar£„ 

Consider  the  element  configuration  shown  in  Figure  1.  A  check  on 
the  calculated  displaceraents  can  be  made  hy  imposing  the  continuum  derived 
boundar}/  conditions  on  a  LHP.  region  where  the  LHPv  is  used  to  simulate  a  com¬ 
pletely  homogeneous  region.  Displacements  at  the  node  points  were  calculated 
and  used  to  calculate  average  stresses  within,  the  elements.  The  computed 
values  of  stress  and  displacement  were  then  compared  with  values  calculated  by 
the  continuum,  solution  for  the  interior  of  the  region.  Agreement  betieeen  the 
two  solutions  \'jas  found  to  be  quite  good.  The  check  thus  provided  one  necessary 
verification  of  both  the  LKB.  element  itself  and  of  the  finite  element  assemblage 
and  solution  procedure  as  well. 

xAs  described  in  the  above ^  all  of  the  energy  released  v/ithin  the 
elem-ent  at  the  crack  tip 3  if  an  increm.ental  am.ount  of  crack  extension  were  to 
occur.,  should  equal  the  strain  energy-release  rate 3  G.  To  check  thlSg  cal- 
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culations  were  made  using  the  exact  near  crack  tip  displacement  field  for  an 
elastic  isotropic  material.  The  energy-release  rate  obtained  from  this  calcula¬ 
tion  is  denoted  by  G*  to  distinguish  it  from  the  exact  theoretical  value  G.  The 
necessity  for  introducing  this  distinction  is  as  follows. 

The  quantity  G  can  be  formally  defined  in  terms  of  a  virtual  crack 
extension  as  the  total  change  in  energy  of  the  body  per  unit  area  of  crack 
extension;  cf,  Equation  (13).  The  quantity  G*,  on  the  other  hand,  reflects 
only  the  change  in  energy  of  the  element  near  the  crack  tip.  Consequently,  G* 
does  not  include  the  change  in  energy  of  the  remainder  of  the  body  as  the  crack 
extends.  In  order  to  account  for  the  contribution  arising  from  the  change  in 
energy  of  the  remainder  of  the  structure,  a  term  denoted  by  6 G *  was  formulated. 
The  sum  G  *  +  6G*  is  then  taken  as  the  appropriate  approximation  to  G  in  this 
work. 

A  numerical  experiment  was  performed  to  calculate  G*  and  6G*  as  a 
function  of  the  aspect  ratio  Ay/ Ax  of  the  LHR  elements.  The  results  are  shown 
in  Table  II.  It  is  quite  evident  that  the  sum  G*  +  6G*  provides  an  accurate 
approximation  to  G  for  the  elongated  aspect  ratios  that  are  convenient  to  use 
in  the  calculations  on  fiber  composite  materials. 

There  are  two  further  points  of  interest  about  the  results  shown  in 

Table  I,  First,  while  it  may  be  surprising  that  6G*  and  G*  are  functions  only 
of  the  ratio  of  Ay  to  Ax,  this  is  a  direct  consequence  of  the  use  of  rigid 
boundary  conditions  which  put  the  near  tip  displacements  into  the  energy  rela¬ 
tionship  used  to  calculate  G.  The  resulting  expression  contains  only  the 
aspect  ratio,  Ay/Ax.  Second,  while  it  might  be  expected  that  SG*  should  con¬ 
verge  to  zero  as  the  grid  spacing  is  collapsed  (i.e.,  as  both  Ay^O  and  Ax^O), 
it  can  be  seen  that  this  was  not  the  case.  Instead,  convergence  is  obtained 
only  as  the  ratio  of  Ay  to  Ax  becomes  large.  Then,  6G*  approaches  a  value  equal 
to  11  percent  of  G  while  G*  converges  to  value  equal  to  86  percent  of  G. 

A  number  of  different  computations  for  heterogeneous  materials  have 
also  been  made.  These  runs  have  been  made  with  a  model  similar  to  the  one  shown 
in  Figure  3.  Preliminary  results  indicate  that  most  of  the  initially  envisioned 
modes  of  failure  are  exhibited  by  the  model.  A  discussion  of  these  results  is 
given  in  the  next  section  of  this  report. 
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TABLE  I.  AN  EXAMPLE  CALCULATION  OF  THE  STRAIN  ENERGY 

RELEASE  RATE  AS  A  FUNCTION  OF  THE  GRID  SIZE  RATIO 
FOR  A  HOMOGENEOUS  LINEAR  ELASTIC  MATERIAL 


la 

Ax 

G* 

G 

6G* 

G 

(G*+5G*) 

G 

1 

0.71 

0.75 

1.46 

2 

0.72 

0.42 

1.14 

4 

0.77 

0.25 

1.02 

8 

0.81 

0.18 

0.99 

16 

0.83 

0.15 

0.98 

32 

0.84 

0.13 

0.97 

64 

0.85 

0.12 

0.97 

128 

0.86 

0.11 

0.97 

256 

0.86 

0.11 

0.97 
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EXAMPLE  COMPUTATIONAL  RESULTS  FOR 
FRACTURE  OF  FIBER- REINFORCED  COMPOSITES 


A  number  of  computations  have  been  performed  using  the  analysis 

technique  described  in  the  preceding  section  of  this  report.  As  discussed 

in  this  section,  the  results  demonstrate  the  ability  of  the  model  to  exhibit 

most  of  the  actual  failure  mechanisms  in  fiber-reinforced  composite  materials. 

These  include  fiber— matrix  debonding,  fiber  bridging,  matrix  bridging,  and 

matrix  crazing.  The  only  important  micromechanical  mechanism  that  the  model 

* 

currently  does  not  explicitly  represent  is  fiber  pull-out. 

Calculations  With  Arbitrarily  Varied  Rupture  Properties 

The  following  describes  a  series  of  example  calculations  in  which 
the  various  fracture  mechanisms  are  made  to  manifest  themselves.  The  re¬ 
sults  shown  here  were  obtained  by  varying  each  of  the  constituent’s  rupture 
properties  and,  in  some  instances,  the  constituent  elastic  properties,  over  a 
wide  range  of  values.  Hence,  the  actual  values  selected  to  perform  these 
example  computations  are  not  totally  realistic.  The  intent  here  is  to  provide 
a  qualitative  verification  that  the  model  is  capable  of  achieving  its  primary 
function  rather  than  to  make  realistic  predictions. 

Unless  otherwise  stated,  the  elastic  constants  used  in  the  calcula¬ 
tions  are  those  given  in  Table  il. 

TABLE  Ill.  ELASTIC  PROPERTIES  USED  IN  THE  SIMULATION 
OF  A  GRAPHITE  EPOXY  COMPOSITE 


Constituent 

E 

Elastic 

Modulus  (ksi) 

V 

Poisson’s 

Ratio 

Fiber 

28,000 

0.3 

Matrix 

495 

0.3 

Interface 

495 

0.3 

*  As  discussed  in  the  Recommended  Future  Research  section  of  this  report, 
steps  to  incorporate  both  fiber  pull-out  and  interply  delamination  can 
be  performed  in  a  straightforward  manner  in  subsequent  work. 
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These  properties  were  intended  to  be  nominally  equivalent  to  those  of  a 
graphite  epoxy  composite  with  a  volume  f rad  ion  of  fiber  equal  to  76  percent. 
Consistent  with  these  values  are  the  following  properties  of  the  bulk  composites 


=  18,000  ksi,  -  0.25 

E  =  690  ksi,  G  -  3,000,000. 
2 


Note  that  the  constituents  were  assum.ed  to  be  isotropic.  However,  the  model 
does  not  require  this. 

The  first  computation  to  be  discussed  is  shovm  in  Figure  6.  Tnis 
result  illustrates  the  fiber-matrix  debonding  mechanism,  a  mechanism  that  occurs 
in  a  large  number  of  cases.  Typically,  a  crack  will  advance  through  the  matrix 
and  then  vertically  separate  or  axially  split  the  interface.  In  this  exaurple , 
the  interfacial  elements  were  made  relatively  weak  while  the  fiber  and  matrix 
material  were  given  relatively  average  strengths.  The  loading  in  this  exam;plc 
was  purely  Mode  1.'  It  was  found  that  axial  splitting  is  even  more  prevalent 
when  Mode  II  loads  are  applied. 

Matrix  crazing  typically  occurs  when  a  crack  advancing  through  the 
matrix  reaches  a  very  stroiig  stiff  fiber.  Figures  7  and  8  sliow  typical  examples 
If  the  fiber  does  not  break,  a  number  of  local  events  occur  In  the  matrix  thfit 
seem  to  approximate  matrix  crazing.  Typically,  the  crazing  does  not  occur  on 
the  uncracked  side  of  the  fiber.  The  reason  is  that  tiie  displacements  on  the 
cracked  side  of  the  fiber  are  "locked  in"  by  the  adjacent  highly  stretched 
fiber.  Consequently,  crazing  usually  occurs  v/hen  stroiig  very  stiff  fiber 
properties  are  inserted  into  the  model  together  with  average  strength  and  stiff¬ 
ness  properties  for  the  matrix  and  interface. 

Matrix  bridging  is  shown  in  Figure  9.  This  phoriornenou  typically 
occurs  when  the  fibers  are  considered  to  be  stiff,  but  v;eak.  lu  these  in¬ 
stances  the  matrix  and  interface  are  capable  of  withstanding  biglior  elonga¬ 
tion  than  the  fiber.  So,  they  will  remain  intact  while  tlie  fiber  cracks.  in 
the  example  shown  in  Figure  9,  some  additional  in ter facial  response  can  also 
be  noticed. 

Unless  otherwise  stated,  in  the  calculations  described  in  this  section  of  the 
report,  the  applied  loading  was  a  tensile  normal  stress  in  the  direclion 
parallel  to  the  fibers. 
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EXAMPLE  CALCULATION  WITH  WEAK  INTERFACE  SHOWING  MATRIX-FIBER  DEBONDING 
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FIGURE  7,  EXAMPLE  CALCULATION  WITH  STRONG  STIFF  FIBERS  AND  WEAK  INTERFACE  SHOWING  MATRIX  CRAZING 
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FIGURE  8.  EXAMPLE  CALCULATION  WITH  STRONG  STIFF  FIBERS  SHOWING  MATRIX  CRAZING 
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Fiber  bridging  is  more  prevalent  when  the  fiber  modulus  does  not 
greatly  exceed  the  matrix  modulus.  (In  graphite  epoxy,  the  modulus  ratio  is 
on  the  order  of  50  and  fiber  bridging  does  not  then  seem  to  occur.)  Figure  10 
demonstrates  this  behavior.  Fiber  bridging  can  also  be  related  to  the  micro- 
structural  geometry.  That  is,  in  computations  performed  with  relatively  thin 
fibers,  this  mechanism  becomes  more  dominate. 

It  should  be  emphasized  that  the  examples  given  in  the  above  deliber¬ 
ately  used  exactly  the  same  LHR  configuration  and  applied  loading.  The  program 
has  the  flexibility  to  consider  different  loading  conditions  and  LHR  geometries. 
However,  this  option  was  not  exercised  here  because  of  the  complications  that 
would  be  added  to  the  interpretation  of  the  results.  That  is,  it  would  be 
difficult  to  separate  the  effects  due  to  the  LHR  geometry  and  loading  from 
those  due  to  variations  in  the  rupture  strength  and  modulus.  Such  calculations 
will  be  deferred  until  further  refinements  have  been  incorporated  into  the 
model.  Again,  the  purpose  of  the  computations  given  in  Figures  6-10  is  to  give 
a  qualitative  demonstration  that  the  model  is  capable  of  coping  with  the  micro¬ 
mechanical  failure  processes  involved  in  the  fracture  of  composite  materials, 
not  to  produce  precise  quantitative  results. 

Finally,  note  that  the  relative  load  levels  at  which  each  individual 
fracture  event  occurs  is  recorded.  (These  levels  are  relative  to  the  load 
level  at  which  the  initial  rupture  event  occurs.)  It  can  be  seen  that  the 
load  level  must  ordinarily  be  increased  to  obtain  additional  rupture  events, 
or,  in  other  words,  in  order  to  enlarge  and  propagate  the  crack-tip  damage 
zone.  This  can  be  contrasted  with  the  behavior  of  completely  homogeneous, 
perfectly  brittle  materials  represented  by  LEFM  where  crack  extension,  once 
initiated,  would  continue  iri  a  catastrophic  manner  under  a  constant  load  level. 
The  initial  stages  of  the  fracture  event  in  fiber  composite  materials  therefore 
can  obviously  be  characterized  in  terms  of  a  stable  growth  process.  This  sug¬ 
gests  that  a  modification  of  the  crack  growth  resistance  curve  (R  curve)  ap¬ 
proach  developed  for  the  fracture  of  ductile  materials  could  be  useful  for 
studying  fractures  in  fiber  composite  materials.  Such  an  approach,  it  might 
be  mentioned,  has  already  been  suggested  in  the  literature;  for  example,  by 
Gaggar  and  Broutman  [9].  However,  these  approaches  are  usually  semiempirical 
in  nature,  relying  heavily  on  LEFM  concepts.  The  present  approach,  in  con¬ 
trast,  should  be  able  to  make  a  direct  prediction  of  the  crack  growth  resis¬ 
tance  parameter  purely  from  fundamental-level  considerations. 
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FIGURE  10.  EXAMPLE  CALCULATION  WITH  STIFF  MATRIX  SHOWING  FIBER  BRIDGING 
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Comparison  of  Calculated  Results 
with  VPISU  Experimental  Data 

The  computations  described  in  the  preceding  section  demonstrate 
that  the  model  developed  in  this  report  displays  the  various  micromechanical 
failure  mechanisms  actually  exhibited  by  fiber  reinforced  composites.  A 
more  quantitative  verification  is  also  possible.  This  can  be  done  by  com¬ 
paring  the  predictions  of  the  model  with  the  experimental  results  obtained 
in  the  concurrent  NASA-Ames  sponsored  research  at  Virginia  Polytechnic 
Institute  and  State  University  (VPISU)  by  Brinson  and  Yeow  [10].  They  have 
measured  the  strengths  of  both  unidirectional  and  angle  ply  graphite/epoxy 
composites  using  unnotched,  single  edge  notched,  and  double  edge  notched  speci¬ 
mens  with  the  crack  introduced  at  various  angles  to  the  fiber  direction.  In 
this  section,  the  model  will  be  used  to  estimate  the  strengths  of  some  of  the 
Brinson-Yeow  unidirectional  notched  tests  using  input  values  inferred  from 
their  unnotched  tests. 

A  precise  prediction  of  failure  loads  is  not  to  be  expected  at  the 
present  stage  of  this  research  for  several  reasons.  First,  because  of  the 
stable  damage  growth  that  precedes  fracture  in  composites,  computations  per¬ 
formed  using  rigid  boundary  conditions  (which  artifically  constrain  the  damage 
growth  process)  will  not  be  realistic.  The  predictions  of  the  initial  local 
failure  event — the  threshold  of  stable  damage  growth — should  be  reasonably 
well  predicted,  however.  This  load  level  will  therefore  be  taken  as  the 
prediction  of  the  model  to  be  compared  with  the  experimental  measurements. 

It  should  be  recognized  that,  because  of  the  neglect  of  the  stable  growth 
regime,  these  should  underestimate  the  actual  failure  loads. 

A  second  reason  for  the  lack  of  precision  inherent  in  the  present 
predictive  capability  is  that  the  constituent  properties  needed  in  the  model 
have  not  been  properly  determined  as  yet.  As  in  the  preceding  section,  linear 
elastic-perf ectly  brittle  behavior  can  be  assumed  with  handbook  values  being 
taken  for  the  elastic  properties.  Rupture  properties  are  not  as  readily 
available,  however.  To  circumvent  this  difficulty,  the  experimental  results 
on  unnotched  specimens  obtained  by  Brinson  and  Yeow  can  be  used.  Their  data 
on  the  strengths  of  unnotched  coupons  pulled  to  failure  with  the  load  in  the 
fiber  direction  (0  =  0°)  and  normal  to  the  fiber  direction  (6  =  90°)  are  given 
in  Table  III.  Relying  on  their  observation  that  matrix  failure  occurred  in 
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TABLE  III.  CRITICAL  VALUES  FOR  MATRIX  FAILURE 
IN  GRAPHITE  EPOXY  COMPOSITES 


Angle 

Between 

(a) 

Experimental  Values 

Critical 

Load  and 

Fiber 

Directions 

Elastic 

Modulus 

(ksi) 

Fracture 

Stress 

(ksi) 

Fracture 

Strain 

(%) 

Strain  Energy 
Density 
(  in .  lb  /  in .  3 ) 

0° 

18,200 

154.7 

0.82 

576 

90° 

1,700 

6.1 

0.44 

13.6 

.a)  Data  of  Brinson  and  Yeow 
a  loading  rate  of  2  x  10“ 

[10]  on  unidirectional 
inches /min. 

unnotched 

specimens  for 
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virtually  all  cases,  the  critical  strain  energy  density  quantity  can  be  cal- 
1  2 

culated  (i.e.,  — Ec  )  for  use  in  the  model.  These  figures  are  given  in  the 
last  column  in  Table  III. 

As  discussed  on  page  11  of  this  report,  the  applied  loads  acting 
on  the  body  are  communicated  to  the  crack  tip  elements  via  the  anisotropic 
elastic  stress  intensity  factors.  For  test  specimens,  such  as  the  single  and 
double  edge  notched  configurations  used  by  Brinson  and  Yeow,  there  is  a  term 
in  the  stress  intensity  factor  that  contains  the  effects  of  the  finite  geometry. 
This  term  is  not  the  same  for  anisotropic  and  isotropic  bodies.  However,  in 
view  of  the  approximate  nature  of  the  present  calculations,  the  refinement 
added  through  the  use  of  the  rather  complicated  anisotropic  expression  was 
not  believed  to  be  warranted.  Hence,  the  isotropic  expressions  were  used. 

For  Mode  I  loading,  these  have  the  form. 


=  o  Y  (2^) 

where  a  is  the  applied  stress  normal  to  the  crack  plane,  a  is  the  crack  length 
and  Y  is  a  dimensionless  function  of  the  specimen  geometry.  For  the  double 
edge  notched  configuration 


Y  =  1.99  +  0.76  (-J-)- 8.48  (|-)  +  27.36  (4-)  ,  (25) 

while  for  the  single  edge  notched  configuration 

Y  =  1.99-0.4l(^)  +  18.7  (^)'-  38.48  +  53.85  {-^)  ,  (26) 

where  W  is  the  plate  width. 

For  a  given  crack  length  and  load-fiber  orientation  in  either  the 
single  or  double  edge  notched  configuration,  the  load  corresponding  to  the 
threshold  of  damage  is  determined  by  a  single  calculation.  This  is  due  to 
the  linear  elastic  material  behavior  assumed  in  the  current  model.  That  is, 
a  computation  can  be  performed  for  an  estimated  value  of  K^.  The  solution 
can  then  be  scanned  to  determine  the  highest  strain  energy  density  in  any 
matrix  element.  By  '‘scaling  up"  the  solution  to  force  this  value  to  match 
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the  critical  value  given  in  Table  IV,  the  value  of  corresponding  to  the 
threshold  of  damage  is  determined.  The  final  step  is  to  use  Equation  (24) 
to  calculate  the  applied  stress  for  direct  comparison  with  the  experimental 
results . 

Comparisons  of  the  calculated  results  as  a  function  of  crack  length 
with  the  Brinson-Yeow  results  on  unidirectional  composites  are  shown  in 
Figures  11  and  12  for  the  single  and  double  edge  notch  configurations,  re¬ 
spectively.  Note  that  a  semilog  format  is  used  (so  that  both  the  0°  and  the 
90°  fiber-load  angle  results  can  be  shown  on  the  same  plot)  and  this  may  make 
the  agreement  appear  to  be  better  than  it  really  is.  Nevertheless,  in  view  of 
the  remarks  made  in  the  preceding  section,  the  prediction  is  reasonably 
accurate  and,  as  expected,  generally  provides  a  lower  bound  to  the  experimental 

data.  Note  also  that  the  unnotched  result  (-—  =  0)  is  included  on  the  theo- 

w 

retical  curve  to  emphasize  that  this  point  was  used  to  construct  the  curve. 

Finally,  comparisons  with  experimental  data  can  also  be  made  for  the 
case  of  cracks  introduced  at  an  angle  to  the  fiber  direction.  An  example 
computation  is  shown  in  Figure  13.  The  threshold-of-damage  load  level  calcu¬ 
lations  are  made  as  above  except  that  both  the  Mode  I  and  Mode  II  stress 
intensity  factors  now  are  involved  in  the  computation.  These  have  been  taken 
as 


i  2 

=  oCL^  Y  sin  ^ 

and 


1 

Kii  =  oa^'  Y  cos  6  sin  (j)  ,  (27) 


where  (f)  is  the  angle  between  the  crack  plane  and  the  fiber  direction.  A 
comparison  between  the  predicted  results  and  the  Brinson-Yeow  experimental 
results  as  a  function  of  (|)  for  a  0°  load-fiber  angle  is  shown  in  Figure  14. 
Again,  it  can  be  seen  that  the  model  provides  a  reasonable  lower  bound  to 
the  experimental  results. 
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FIGURE  11.  COMPARISON  OF  PREDICTED  STABLE  GROWTH  THRESHOLD  WITH 

EXPERIMENTAL  FAILURE  LOADS  FOR  SINGLE  EDGE  NOTCH  SPECIMENS 


Remote  Applied  Stress ,  ksi 


FIGURE  13.  EXAMPLE  CALCULATION  FOR  A  FIBER  COMPOSITE  WITH  AN  ANGLE  CRACK 
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RECOMMENDED  FURTHER  RESEARCH 


Being  two  dimensional,  the  model  described  in  this  report  is  so  far 
limited  in  application  to  unidirectional  laminates  and  can  only  reflect  the 
local  damage  modes  that  occur  in  plane  deformation.  In  subsequent  work,  it 
would  be  desirable  to  extend  the  model  from  two  dimensions  to  three  dimensions 
to  treat  angle  ply  laminates.  In  addition  to  treating  the  effect  of  different 
layups,  new  modes  of  damage  (e.g.,  interply  delamination,  fiber  pull-out) 
should  also  be  included. 

A  straighf orward  extension  of  the  model  could  be  made  by  devising 
a  three-dimensional  LHR,  But,  the  number  of  storage  locations  required  for 
the  three-dimensional  grid  would  likely  exceed  the  fast  access  storage  of  most 
computers.  To  circumvent  this  situation,  the  present  two-dimensional  model 
could  be  extended  to  angle  ply  laminates  by  considering  laminates  with  an  LHR 
only  in  a  single  ply.  Figure  15  shows  this  concept.  In  the  simplest  case, 
the  LHR  can  be  confined  to  the  most  critical  ply  as  the  load  level  is  increased. 
In  general,  there  can  be  an  LHR  in  each  ply  with  attention  shifted  from  one  ply 
to  another  in  turn.  In  either  case,  it  would  be  necessary  to  assume  that  the 
interply  interactions  are  not  affected  by  the  precise  details  of  the  local 
damage  occurring  in  each  ply. 

In  Figure  15a,  a  laminate  containing  a  flaw  is  shown  with  generalized 
loading  conditions.  In  this  the  first  level  of  modeling,  all  of  the  plys  are 
considered  as  homogeneous  orthotropic  layers.  Displacements  can  be  calculated 
in  this  model  for  any  load  level  and  applied  to  the  boundaries  of  the  local 
homogeneous  region  in  each  ply.  Then,  just  as  in  the  two-dimensional  model, 
the  LHR  shown  in  Figure  15b  can  be  monitored  to  determine  the  local  rupture 
events  and  the  extent  of  the  damage  zone  at  that  load  level.  When  significant 
crack  growth  has  occurred  in  one  or  more  plys,  the  first  level  model  must,  of 
course,  be  made  to  reflect  this.  But,  the  basic  approach  would  not  be  changed 
thereby. 

Two  distinct  techniques  must  be  developed  in  order  to  obtain  an 
accurate  set  of  boundary  displacements  for  each  LHR  when  extensive  damage 
occurs  in  that  or  in  neighboring  plys.  The  first  is  to  incorporate  the  effect 
of  the  stiffness  change  in  the  LHR  into  the  first  level  model  of  the  laminate. 

loading  conditions  would  then  be  obtained  from  this  modified  first  level 
model  (i.e.,  one  with  locally  reduced  stiffnesses)  to  give  the  new  boundary 
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conditions  for  the  LHR.  In  the  second  technique,  the  influence  of  the  stiff¬ 
ness  of  the  laminate  away  from  the  crack  tip  will  be  directly  transferred  to 
the  LHR  through  what  can  be  called  a  "flexible  boundary  condition"  technique. 

This  approach  can  be  likened  to  placing  a  series  of  springs  on  the  LHR  boundary 
with  stiffnesses  taken  from  the  first  level  model.  The  actual  technique  is 
more  sophisticated  than  this,  however,  as  follows. 

In  the  work  described  so  far,  the  peripheral  elements  in  the  LHR 
have  been  positioned  according  to  the  displacement  field  given  by  a  completely 
linear  elastic  continuum  solution  for  the  given  crack  tip  location  and  applied 
loads.  This  is  not  strictly  correct  even  when  the  displacements  are  periodically 
updated  as  the  crack  extends  in  the  LHR.  Even  in  the  absence  of  local  damage, 
the  highly  nonlinear  inhomogeneous  nature  of  the  crack-tip  region  in  a  com¬ 
posite  will  cause  significant  departures  from  the  continuum  displacements.  A 
scheme  for  estimating  these  displacements  based  on  a  determination  of  the  equi¬ 
librium  state  between  the  LHR  and  the  continuum  can  be  obtained  by  adapting 
the  technique  reported  in  Reference  3.  Such  a  scheme  has  been  referred 
to  as  the  "flexible  boundary  condition"  approach.  The  approach  used  so  far  can 
be  called  the  "rigid  boundary  condition"  approach  because,  even  with  periodic 
updating  to  reflect  the  progress  of  the  crack  through  the  LHR,  the  LHR  does 
not  directly  affect  the  continuum  region  surrounding  it.  The  flexible  boundary 
condition  approach,  in  contrast,  accounts  for  the  interaction. 

Finally,  the  range  of  modes  of  damage  should  be  extended  over  those 
of  the  two-dimensional  model  to  include  interply  delamination  and  a  "free- 
edge"  effect.  The  general  approach  could  be  similar  to  that  described  in 
References  10  and  11.  Attention  should  first  be  focused  on  through- the- thickness 
cracks  in  laminated  plates  under  tension  with  part-through  flaws  being  considered 
later.  Of  most  importance,  the  three-dimensional  model  of  an  angle  ply  laminate 
should  permit  arbitrarily  varied  multidirectional  layups  to  be  explicitly  con¬ 
sidered.  The  properties  for  the  ply  moduli  could  be  obtained  from  laboratory 
investigations,  e.g.,  in  the  program  being  conducted  for  NASA- Ames  at  Virginia 
Polytechnic  Institute  and  State  University. 
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